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Abstract 

The spreading of liquid nanodroplets of different initial radii Rq is studied using molecular dy- 
namics simulation. Results for two distinct systems, Pb on Cu(lll), which is non-wetting, and a 
coarse grained polymer model, which wets the surface, are presented for Pb droplets ranging in size 
from ~ 55 000 to 220 000 atoms and polymer droplets ranging in size from ~ 200 000 to 780 000 
monomers. In both cases, a precursor foot precedes the spreading of the main droplet. This precur- 
sor foot spreads as r 2 j(t) = 2D e fft with an effective diffusion constant that exhibits a droplet size 
dependence -D e // ~ ' ■ The radius of the main droplet r^t) ~ R^ 5 in agreement with kinetic 
models for the cylindrical geometry studied. 
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Numerous practical problems are determined by the wetting of surfaces by liquids includ- 
ing adhesion, lubrication, coatings, plant protection, and oil recovery. The spreading of a 
droplet on a surface proceeds as the droplet balances the interfacial tensions between the 
solid, liquid, and vapor phases. Completely wetting droplets are known to form a monolayer 
film or precursor foot on the surface that spreads ahead of the droplet. For "high-energy" 
surfaces, spreading includes the formation of molecule-sized terraces along with the precursor 

foot QQ. 

Although the spreading behavior of droplets on surfaces has been extensively studied, 
droplet size effects on the spreading dynamics are often ignored. These size effects can have 
a large influence on the wetting behavior, particularly when using nanoscale droplets such as 
in MEMS and microfluidic devices. For these devices, droplet size can strongly affect both 
the manufacturing speed, such as in microcontact printing, and the device performance. 

For a spreading liquid droplet, the energy dissipation mechanics have been classified by 
de Gennes. The total dissipation in the spreading droplet can be expressed as a sum of 
three distinct dissipation mechanisms, TS W + TT1/ + TS/ [3]. In this equation, TH W is the 
contribution due to viscous dissipation in the bulk of the droplet, TE f is the contribution due 
to viscous dissipation in the precursor foot, and TE; is the contribution due to adsorption of 
liquid molecules to the surface at the contact line. The adsorption mechanism is expected 
to dominate for low viscosity systems at short times while the bulk viscous dissipation 
mechanism takes over at later times The dissipation from the foot is less understood, 
though in the simulations presented here we will show that this dissipation does not play a 
role in the spreading of the droplet. 

The size dependence of the bulk droplet spreading rate is present in models of droplet 
spreading. In the molecular kinetic theory of liquids 0,0], the energy dissipation occurs at 
the contact line. In the linearized version of this model, the contact radius r b (t) of the bulk 
region of the droplet scales with the droplet volume according to r&(i) ~ Rq at late times for 
spherical droplets [3| and 77, (t) ~ Kq 5 for cylindrical droplets The hydrodynamic model 
[3, is based upon the solution of the equations of motion and continuity for the droplet 
and assumes that energy dissipation occurs via viscous dissipation in the bulk. For the 
hydrodynamic model, the contact radius of the bulk region scales according to r b (t) ~ Rq^ 10 
at late times for spherical droplets Q and r&(t) ~ R^ 7 for cylindrical droplets jfj. de Ruijter 
et al. 0| derived a combined model by including the energy dissipation mechanisms from 
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both the linearized kinetic and hydrodynamic models. This model does not consider the 
viscous dissipation in the precursor foot. For a cylindrical droplet, the combined spreading 
model is 



dr b {t) 7 (6 



dt Co , 2y(r b (t)-a) sin 2 e I s j n Q s j n Q 

2 ~ r 6 (t)(0-sin0cos0) v 



(1) 



d6 _ (6 -smOcosO\ 1/2 ( sin 3 6 \ 1 dr b (t) 

dt~{ a J \;° e- sine cose J ~HT' (2) 

where the cross-sectional area A ~ R^, 7 is the liquid/vapor surface tension, ( is the 
friction coefficient, r\ is the bulk fluid viscosity, a is the hydrodynamic cutoff, and #0 is the 
equilibrium contact angle. The linearized kinetic and hydrodynamic models are obtained 
from Eq. [T] by setting rj and (q to zero, respectively. Although this produces the linearized 
version of Blake's kinetic model, experiments 0, U, 3 and simulation 0, Q, have 



shown that it works quite well. These models are written for droplets with a finite 9 , but 
they can also be applied to wetting droplets by setting = 1 in Eq. [TJ 

Droplets often spread by extending a precursor foot ahead of the main droplet. This has 



been observed in experiments 



id 13| as well as simulation 0, [3, Q| ■ This foot grows 



diffusively, though to the best of our knowledge there are no predictions for the dependence 
on the droplet size. Droplet spreading experiments have observed "terraced spreading" where 
multiple layers spread on top of the precursor foot. In simulations of droplet spreading, the 
precursor foot is present, but terraced spreading has not been observed. This is probably 
due to the relatively short duration of the simulation runs. 

Although droplet size effects are often ignored when studying spreading dynamics, we 
show that the spreading rate of both the precursor foot and the bulk change with droplet 
size. Here, we present MD simulations for two very different systems, a coarse-grained model 
of polymer nanodroplets in the wetting regime and an explicit atom model of Pb on Cu(lll), 
which is non- wetting, to study the dependence of the spreading rate of both the droplet and 
the precursor foot on the initial radius i?o of the droplet. Our results demonstrate that the 
observed behavior is a general phenomenon and does not depend on the system specifics. In 
both cases, the vapor pressure is low so that spreading does not occur via vaporization and 
condensation. 

For polymer chains, the polymer is represented by spherical beads of mass m attached 
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by springs, which interact with a truncated Lennard- Jones (LJ) potential, 



U LJ (r) 



is 



(?) 
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r < r r 

(3) 
r > r r 



where e and a are the LJ units of energy and length and the cutoff r c = 2.5a. The monomer- 
monomer interaction e is used as the reference and all monomers have the same diameter a. 
For bonded monomers, we apply an additional potential where each bond is described by 
the FENE potential |20j with k = 30 e/a 2 and Rq = 1.5 a. The substrate is modeled as a 
flat surface since it was found previously (3| that with the proper choice of thermostat, the 
simulations using a flat surface exhibit the same behavior as a realistic atomic substrate with 
greater computational efficiency. The interactions between the surface and the monomers in 
the droplet at a distance z from the surface are modeled using an integrated L J potential with 
the cutoff set to z c = 2.2a Extending the range of this surface interaction to infinity 

increases the spreading rate of the precursor foot and slightly increases the spreading rate 
of the bulk. Aside from shifting the wetting transition to a lower energy, the qualitative 
spreading behavior is identical to the z c = 2.2a case. Specifically, the absence of terraced 
wetting is probably not due to the finite interaction range of the surface used in simulations 
but to the length of the simulations. Here we present results for e w = 2.0 e which for N = 10 
is above the wetting transition e c w = 1.75 e. 

We apply the Langevin thermostat to provide a realistic representation of the transfer of 
energy in the polymer droplet. The Langevin thermostat simulates a heat bath by adding 
Gaussian white noise and friction terms to the equation of motion, 

m i -^ = -AU i -ma L -^- + W i (t), (4) 

where rrii is the mass of monomer i, 7^ is the friction parameter for the Langevin thermostat, 
— AUi is the force acting on monomer i due to the potentials defined above, and Wj(t) is a 
Gaussian white noise term. Coupling all of the monomers to the Langevin thermostat would 
have the unphysical effect of screening the hydrodynamic interactions in the droplet and not 
damping the monomers near the surface stronger than those in the bulk. To overcome this, 
we use a Langevin coupling term with a damping rate that decreases exponentially away 
from the substrate |2l]. We choose the form 71,(2) = 7£exp (a — z) where Yl i s the surface 
Langevin coupling and z is the distance from the substrate. In this paper, we present results 



for Yl = 3-0 and 10.0 r _1 . The larger Yl corresponds to an atomistic substrate with large 
corrugation and hence large dissipation and slower diffusion near the substrate. 

For Pb on Cu(lll), interactions are described via embedded atom method (EAM) inter- 
atomic potentials wherein the energy for N atoms is [22 1 

N l 

^ = E[%) + 5 EW ( 5 ) 

i=i z j# 

In Eq. El Pi is the electron density at atom i, pi = J2j^i Pj a ( r ), where pj a {r) is the spherically 
symmetric electron density contributed by atom j, a distance r from i. Fi(pi) is the energy 
associated with embedding atom i into an electron density pi and <f>ij(r) is a pair potential 
between atoms i and j. The many-body nature of Fi(pi) in Eq. El makes the EAM superior to 
pair potentials for describing bonding in metal systems. The interactions for Cu, Pb, and the 
cross-term between them were previously parameterized |2i 



24123- The Cu(lll) substrate 



was described via an explicit atom description with dimension in the surface normal direction 
equal to four times the potential cutoff r c = 5.5 A. The substrate was equilibrated at the 
proper lattice constant prior to joining with the drop. Atoms in the 2r c planes furthest from 
the surface are held rigid and the rest are permitted to relax according to MD equations of 
motion throughout all simulations. Because the substrate is represented atomistically for 
Pb(l) on Cu(lll), we thermostat only atoms in the substrate using a Nose-Hoover thermostat 
algorithm. 

All of the droplets presented here are modeled as hemicylinders as described previously 
The system is periodic in the y direction with length L y and open in the other two 
directions. This allows a larger droplet radius to be studied using the same number of 
monomers than in the spherical geometry. Polymeric droplets have initial droplet radii of 
#0 = 50, 80, and 120 a, cross-sectional areas A = 2690, 7490, and 18 200 a 2 and a total size 
N = 200 000, 350 000, and 780 000 monomers, respectively, with L y = 60 a for R Q = 50 a 
and L y = 40 a for i? > 80 a. Pb(l) droplets are studied with R = 20, 30, and 40 nm, 
A = 630, 1400, and 2500 nm 2 and N = 55 000, 122 000, and 220 000 atoms in the drop, 
respectively, with L y = 27 A. In each case, keeping L y < R or L y fh R supresses any 
Rayleigh instability. The substrates for the three drop sizes in the Pb(l) on Cu(lll) systems 
contained = 85 000, 128 000, and 170 000 Cu atoms. 

The equations of motion are integrated using a velocity- Verlet algorithm. For polymer 
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spreading, we use a time step of At = 0.01 r where r = a ( -j) . The simulations are 



Figure 1: (color online) Cross-section snapshots from the simulations. A droplet composed of 
chain length N = 10 polymers with initial radius Ro = 80 a, thickness L y = 40 a, e w = 2.0 e, and 
Y L = 3.0 r" 1 at t = 71 600 r (top). A Pb(l) droplet with R = 40 nm and L y = 27 A on Cu(lll) 
at t = 3.6ns (bottom). The snapshots show substrates that are 750 a and 230 nm in length for the 
polymer and Pb droplets, respectively. 



performed at a temperature T = e/ks using the LAMMPS code |2J. For Pb(l) on Cu(lll), 
At = 1 fs, T = 700 K. and the code paradyn bjj was used. 

n 

For the polymeric system, it was shown previously [13] by measuring the equilibrium 
contact angle that finite system size effects become negligible for droplets containing 50 000 
monomers or more. Here, we study the size dependence of the spreading behavior for droplets 
substantially larger than this finite size limit by using a minimum droplet size of 200 000 
monomers in order to simultaneously study the bulk and precursor foot regions. This is 
shown in Fig. [2 which shows the foot extending beyond the bulk region. For the monomeric 
liquid system Pb on Cu(lll), computational requirements for modeling the substrate are 
more significant so we use a minimum drop size N ~ 55 000 atoms. In a prior work 
wetting of Pb(l) drops on Cu(100) and Cu(lll) was simulated using this model. For Pb(l) on 
Cu(lll), a precursor film rapidly advanced ahead of the main drop in a diffusive manner (see 
Fig. Q]); furthermore, the drop reached an equilibrium but finite contact an gle = 33° on 
top of the prewetting film, which is in excellent agreement with experiment |28l l29l. l30i. l3l|. 
This system exhibits negligible exchange of atoms between the liquid and solid, permitting 
evaluation of system size effects for non-reactive wetting in the case of a monomeric liquid 
with low vapor pressure. 

For the simulations presented here, we extract the instantaneous contact radius r&(£) and 
contact angle 6{t) every 400 r for the polymer system and every 4ps for the metal system 
according to the procedure described previously P, To demonstrate the influence of the 
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viscous dissipation in the precursor foot on the dynamics of the bulk region, two systems are 
studied where the droplets are placed on substrates prewet with a monolayer of the same 
material to eliminate the simultaneous spreading of the precursor foot. The bulk contact 
radii for both a polymer droplet and a lead droplet are shown in Fig. [21 comparing spreading 
on prewet substrates to spreading on bare surfaces. Adding a monolayer to the substrate 
does not affect the bulk spreading for the polymeric system indicating that the viscous 
dissipation in the precursor foot is negligible for the droplet j^Sfl- Although the Pb on Cu 
system is drastically different than the polymer system, the results are very similar. For the 
lead droplet, the early time behavior shows that the bulk spreading rate is enhanced as the 
foot is formed and begins to extend. Once the foot has extended away from the droplet, the 
bulk contact radius curves become parallel. Plotting the velocities of the precursor foot on 
both the bare and prewet surfaces shows that they are identical after one nanosecond. 

The scaling predictions of the linearized kinetic and hydrodynamic models are applied 
to the bulk contact radius data for e w = 2.0 e and Yl = 3-0 and 10.0 r -1 in Fig. where 
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dividing the bulk contact radius by R for each of three droplet sizes causes the data to 
collapse to a single curve. The R^ 7 scaling of the same data, shown in Fig. Eb, does not fit 
as well because hydrodynamic energy dissipation has only a weak influence on the spreading 
rate for the conditions used in these simulations |8|. Similar results are found for Pb on 
Cu(lll) (see Figs. and b) because here we also have a low viscosity, rapidly spreading 
droplet. 

Previous simulations of polymer droplet spreading Q have shown that the hydrodynamic 
model, Eqs. [U and El with £ = 0, adequately fits the data only for the higher viscosity 
(longer chains) and slower spreading droplets. Here, a consistent improvement of the fit to 
the hydrodynamic model is observed for the three droplet sizes in going from the faster Yl = 
3.0 t -1 to the slower Yl = 10-0 r_1 conditions as well as going from smaller to larger droplets. 
Both the kinetic and combined models fit the data well for all of the droplets alth oug h the 
combined model tends to produce less reasonable values of the fitting parameters 

The spreading of the precursor foot has been measured experimentally by ellipsometry 
and more recently by atomic force microscopy These studies report effective diffusion 
coefficients for the precursor foot without considering the dependence on the droplet size. 
Joanny and de Gennes predicted the height profile of the precursor foot to be proportional 
to l/r 0,0, 

but models relating the precursor foot dynamics to the droplet dimensions, 
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Figure 2: Comparison of bulk droplet spreading rate on prewet surface (solid line), and bare surface 
(dashed line) for (a) N = 10 polymers with e w = 2.0 e and -) S L = 10.0 r _1 and (b) a Pb(l) droplet 
on Cu(lll). 

such as those presented above for the bulk dynamics, are not available. Like the results of 
the bulk contact radius, the size dependence is evident in the precursor foot contact radius 
as well. This is shown in Fig. El where the precursor foot contact radius is divided by the 
initial contact radius raised to the power n where n = 1/2. For the Yl = 10.0 r -1 system, 
the curves in Fig. Eh show the same asymptotic behavior, but differences in the initial 
contact radius cause the curves to be offset by a constant value at later times. For the 
7£ = 3.0 r^ 1 system, Fig. Eh shows the offset is less severe and the curves for the three 
different droplet sizes overlap. For Pb(l) on Cu(lll), Fig. Et> shows results in this system 
are quite similar to what is seen for the Yl = 3.0 r^ 1 case in the polymeric systems. This 
implies that Pb(l) on Cu(lll) corresponds to a lower surface corrugation which, given the 
significant lattice mismatch between Pb and Cu and the high density packing of the (111) 
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Figure 3: Scaling of the bulk contact radius for three droplet sizes based on the initial radii 
Ro = 50 a (solid line), Rq = 80 a (dotted line), and Rq = 120a (dashed line) using the predictions 
of (a) the kinetic model and (b) the hydrodynamic model. Each droplet is composed of chain 
length N = 10 polymers with e w = 2.0 e, 7£ = 3.0 (upper curves) and 10.0 r _1 (lower curves). The 
7£ = 3.0 t -1 curves are shifted upward for clarity. 

surface, is a reasonable result. 

The precursor foot contact radius follows r"j(t) = 2D e fft where -D e // is the effective 

diffusion coefficient. Within the uncertainty of our data, these results are consistent with 

D e ff ~ Rq where x = 0.5 ± 0.05 even though the polymeric systems are completely wetting 

and the metal system is nonwetting. For the Yl = 10.0 r _1 system, D e ff ~ _Rq 65 collapses 

i ii 

the data onto a master curve better than R . This may be due to the fact that the 
high coupling constant reduces the foot spreading rate to that of the bulk droplet, so the 

1/2 

D e ff ~ R ' scaling is not valid when the spreading of the bulk droplet interferes with the 
precursor foot diffusion. At present, there is no theoretical model we are aware of which 
predicts the dependence of -D e // on i? , but it may be because the higher bulk spreading 
rate of the larger droplets pushes the precursor foot outward adding to the driving force of 
the surface interaction. 

In summary, we study the droplet size dependence of nanodroplets spreading in cylindrical 
geometries using molecular dynamics simulation. Our results follow the kinetic model of 
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Figure 4: Scaling of the bulk contact radius of three Pb(l) droplet sizes based on the initial radii 
Ro = 20 ran (solid line), Rq = 30 nm (dotted line), and Rq = AOnm (dashed line) using the 
predictions of (a) the kinetic model and (b) the hydrodynamic model. 

4/5 

droplet spreading, which predicts a R size scaling of the bulk droplet contact radius. The 
bulk spreading rate does not change if the droplet is instead spread on a prewet surface 
consisting of a monolayer of the droplet material. This indicates that, at least in the present 
simulations, the viscous dissipation from the precursor foot is not important for studying 
the kinetics of the droplet. Theories describing the dynamics of the precursor foot do 
not predict a droplet size dependence. The spreading rate of the precursor foot is known 
to be diffusive, but here we show that the effective diffusion coefficient has a droplet size 

1/2 

dependence, -D e // ~ Rq , which can be utilized in the design of surface wetting applications. 
For Pb(l) on Cu(lll), use of a realistic interatomic potential results in non-reactive wetting 
(in agreement with experiment). Furthermore, vapor pressure is low despite this being a 
monomeric liquid. As such, spreading mechanisms for the metal systems are similar to the 
polymeric systems studied and the same size dependence of the spreading rate is found 
despite the fact that these are two very different systems. 

Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin 
Company, for the United States Department of Energy's National Nuclear Security Admin- 
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Figure 5: (a) Polymeric systems precursor foot spreading rate divided by R for Ro = 50 c (solid 
line), Rq = 80 a (dotted line), and Rq = 120 a (dashed line) for e w = 2.0 e and 7£ = 10.0 t -1 
(lower curves) and for 7£ = 3.0 r" 1 (upper curves). Each droplet is composed of chain length 
N = 10 polymers. The 7£ = 3.0 r _1 curves have been shifted upward for clarity, (b) Precursor foot 
spreading rate for Pb(l) on Cu(lll) divided by R^ 2 for R = 20, 30, and 40 ran at T = 700 K. 

istration under Contract No. DE-AC04-94AL85000. 
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